% This function takes a rental rate and spillovers and return the difference between the housing supply
% and housing demand. We search over the rental rate to find the one that sets supply equal to demand
function [share_diff, tax, mean_wage] = find_rents_mto(rs, Ss, Hs, myDist, w_mto)
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha xi xi1 ub lb pi_e H w theta Int_pi_w sigmae elasts...
        phi sigma gamma rts e_curr r_base pe myiter pop_gr vouch_sh part_eqm thresh_flag
    
    % Set variables
    rs = max(rs, 0.0); % The rental rate in A can't be less than B
    mean_wage = squeeze(sum(myDist .* w, 1:2) ./ sum(myDist, 1:2));
    
    wto_share = 0.3; % Share of income going to MTO housing
    tax       = 0.0; % Initial guess for MTO tax

    thresh_grid = [0.0015 0.05 0.10 0.15 0.20 0.25];
 
    if part_eqm == 1
        
        load tax_ss
        thresh_grid = [0.0015 0.05 0.10 0.15 0.20 0.25];
        tax = tax_ss(find(thresh_grid==thresh_flag),myiter);
    end
    
    % Fixed point to find the tax rate that finances the MTO experiment
    for inner_iter = 1:20
        % Calculate choice probabiities
        [mdist, F_es, wprimeS, educ_mat] = choiceProb_ct_par(Ss, rs, tax, w_mto, mean_wage);
        
        % Calculate housing demand
        H_dem = sum(bsxfun(@times, mdist, myDist), 1:3);
        H_dem = squeeze(H_dem);
        H_dem = bsxfun(@times, H_dem, pi_theta);
        H_dem = sum(H_dem, 2);
        
        wMass = sum(bsxfun(@times, mdist, myDist .* w), 1:3);
        wMass = squeeze(wMass);
        wMass = bsxfun(@times, wMass, pi_theta);
        
        % Find neighborhood ranking 
        [~, id_max] = sort(-mean_wage);
        % Calculate MTO voucher cap
        vouch_cap = vouch_sh*rs(id_max(1)); 

        % Tax that clears the below housing market
        tax_11 =  sum(myDist .* (min(rs) - w) .* (w < min(rs)),1:3)/sum(myDist .* w, 1:3);
        % Tax that clears the MTO experiment
        share_A = mdist(:,:, id_max(end), id_max(1), 1)*pi_theta(1)+mdist(:,:,id_max(end), id_max(1), 2)*pi_theta(2);
        tax_12  = sum(share_A .* myDist .* (w <= w_mto) .* min(max(min(rs(id_max(1)) - w*wto_share, vouch_cap ),0), 100*rs(id_max(2))),1:3)/sum(myDist .* w, 1:3);

        tax1 = tax_11*0  + tax_12 ;
         
        if part_eqm == 1
            
            load tax_ss
            thresh_grid = [0.0015 0.05 0.10 0.15 0.20 0.25];
            tax = tax_ss(find(thresh_grid==thresh_flag),myiter);
         end

        err = abs(tax1 - tax);
        if err < 1e-6
            break
        elseif  inner_iter == 20
            disp("Didn't converge in taxes")
            disp(err)
        end 
        tax = tax1;
    end

    % Calculate housing demand
    H_dem = H_dem(:)   ;
    
    % Calculate housing supply
    tmp = housing_supply(rs, Hs);
    
    % Calculate difference accounting for population growth
    share_diff = H_dem.*(pop_gr).^(myiter-1) - tmp;
end